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1.  Introduction 


Our  recent  work  ( 7 ,  2)  in  the  area  of  granular  physics  indicates  that  certain  (ordered  with  respect 
to  size)  arrangements  of  smooth,  metal  spheres  adjacent  to  each  other  in  one-dimensional  (1-D) 
chains  demonstrate  noticeable  shock  absorption.  This  has  been  observed  and  measured 
analytically  (7-5)  and  through  numerical  integration  of  the  equations  of  motion  (7,  6-8). 
Experimental  efforts  (9-12)  have  focused  primarily  on  solitary  wave  propagation  and 
demonstrating  reduced  impact  force  using  tapered  chains  (TCs)  (4).  The  system  is  scalable,  as 
there  does  not  appear  to  be  any  restrictions  on  particle  size.  There  is,  of  course,  limits  imposed 
by  quantum  physics  and  manufacture  of  the  spheres.  Up  to  this  point  in  our  efforts,  only  single 
chains  of  varying  length  have  been  considered.  However,  increasingly  encouraging  results  have 
prompted  us  to  consider  the  collective  disruptive  effects  of  many  chains  as  a  defeat  mechanism 
against  ballistic  shock  and  other  undesirable  transient  pulses. 

The  goal  of  this  report  is  to  summarize  both  published  and  unpublished  results  through 
nonnalized  energy  calculations  and  include  an  analytical  description  of  a  potential  armor  panel. 
The  contents  are  organized  as  follows.  A  small  overview  of  granular  media  and  the  ordered 
systems  of  interest  are  first  discussed.  This  is  followed  by  a  hard-sphere  approximation  and  the 
numerical  solution  to  the  equations  of  motion.  In  addition,  our  mathematical  model  provides 
some  predictive  capability  for  a  certain  range  of  systems.  Past  experimental  arrangements  and 
intended  investigations  are  reviewed.  In  conclusion,  we  introduce  a  potential  annor  panel 
configuration  with  an  estimate  of  the  specific  and  volumetric  absorbed  energy  and  summarize 
the  salient  features  of  the  report. 


2.  Granular  Media  and  TCs 


Granular  media  (13-17) — not  to  be  confused  with  microscopic  grains  in  a  metal — are  a  curious 
thing.  In  general,  its  members  consist  of  discrete  particles  that  can  range  in  size  from 
micrometers  to  meters  and  number  in  quantity  from  several  to  the  uncountable.  The  most  well- 
known  constituents  of  this  group  are  sands  and  powders  and  are  utilized  across  many  disciplines. 
In  this  communication,  components  of  granular  media  may  be  referred  to  as  grains,  particles, 
beads,  or  spheres.  These  media  are  particularly  intriguing  in  that  among  other  things  it  has 
properties  of  liquids  and  solids  and  calls  have  been  made  to  add  it  as  a  fundamental  state  of 
matter  (13).  Within  the  U.S.  Army  Research  Laboratory,  granular  media  have  been  investigated 
for  understanding  wave  propagation  and  its  mathematical  similarity  to  nonlinear  rods  (18). 

While  granular  media  have  been  fundamental  in  saving  lives  for  decades  against  ballistics 
(sandbags),  new  features  are  being  observed  in  the  growing  literature  which  suggest  a  more 
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technological  role  may  be  at  hand.  Particularly  appealing  and  relatively  simple  to  address  is  the 
behavior  of  1-D  systems  where  smooth  elastic  spheres  intersect  at  a  point  in  their  initial  state. 
Even  though  we  may  constrain  the  basis  of  variables  in  this  system  to  just  a  few  parameters,  the 
complexity  of  its  dynamics  is  formidable.  We  call  these  systems  tapered  chains  (TCs). 

We  define  TCs  (figure  1)  as  1-D  granular  arrays  of  elastic  spheres  that  touch  at  a  single  point  in 
their  initial  state  and  grow  to  a  disk  under  compression  in  the  plane  perpendicular  to  the  figure. 
The  chains  can  be  characterized  by  the  number  of  grains,  N,  the  successive  decrease  in  size  of 
the  grains  or  tapering,  q,  and  restitutive  losses,  co.  Restitution  (19-20)  represents  the  expansion 
phase  in  a  collision  where  energy  is  converted  from  elastic  potential  to  kinetic  and  energy  losses 
can  be  incorporated. 


Figure  1.  The  simple  tapered  chain:  N=  10,  q  =  8%,  L  =  70.7  mm,  and  r,=  5  mm. 

It  is  intriguing  and  suggestive  that,  conceptually,  granular  media  can  be  envisioned  as  the  inverse 
of  a  porous  material.  In  a  porous  material,  there  are  gaseous  (air)  voids  within  a  solid  matrix. 

For  dry  granular  media,  solid  voids  (grains)  exist  within  a  gaseous  matrix  (air).  Packing, 
however,  limits  the  analogy;  for  grains,  each  entity  is  in  contact  with  at  least  one  other,  whereas 
voids  are  not  necessarily  in  contact  with  one  another  in  porous  materials.  This  is  an  interesting 
comparison  to  draw  because  the  considerable  amount  of  work  required  to  collapse  the  voids  in  a 
porous  material  makes  it  a  favorable  technology  against  ballistic  shock  (21).  Are  there  then 
opportunities  in  granular  media  for  the  armor  designer  to  exploit? 

The  interaction  potential  between  adjacent  elastic  spheres  was  first  identified  by  Hertz  (22)  circa 
1881.  Modern  and  eloquent  derivations  were  later  performed  by  both  Landau  and  Lifshitz  (23) 
and  Love  (24),  with  a  simplified  order  of  magnitude  approach  by  Leroy  (25).  More  recently, 
Nesterenko  (26)  has  written  extensively  on  the  dynamics  of  granular  chains.  If  one  can  show 
significant  reduction  in  output  energy  while  ensuring  that  the  structure  of  granular  columns  can 
be  maintained,  it  is  quite  possible  to  use  these  TCs  against  propagating  shocks.  For  now,  we 
assume  that  one  remains  within  the  elastic  regime.  Future  work  will  begin  to  relax  this 
restriction. 
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3.  Hard-Sphere  Approximation 


These  approximations  differ  from  the  numerically  simulated  systems  in  two  major  ways.  The 
first  is  that  the  chain  is  not  bounded  by  fixed  rigid  walls.  As  a  result,  energy  will  not  continue  to 
be  transmitted  up  and  down  the  chain.  The  second  is  that  the  potential  becomes  infinite  and  as  a 
consequence,  the  energy  packet  is  only  1  grain  in  width.  The  system  therefore  propagates  energy 
as  independent  collisions.  This  is  congruent  to  the  independent  collision  model  proposed  by  Wu 
(■ 4 ).  In  the  numerical  results,  width  of  the  energy  pulse  is  a  function  of  the  tapering — and  when 
all  grains  are  the  same  size  it  is  about  5  grains  wide.  By  generating  an  iterative  form  of  the 
conservation  equations,  one  can  arrive  at  an  expression  for  the  normalized  kinetic  energy  (KE), 
KEn  =  KEout/KEin .  This  ratio  will  be  the  primary  variable  determining  the  absorptive  quality  of 

TCs. 


3.1  The  Simple  Tapered  Chain 

The  simple  tapered  chain  (STC)  is  displayed  in  figure  1.  To  generate  an  initial  disturbance,  an 
input  velocity,  v,  ,  is  applied  to  the  rightmost  and  largest  grain  with  radius,  rr  It  propagates  to 
the  left,  encountering  an  initially  stationary  grain  of  radius,  rM.  The  radius  of  the  /  +  1  particle 
may  be  reduced  by  q%  from  r  .  This  tapering  q  will  be  constant  along  the  entire  length  of  the 
chain.  During  the  transmission  of  the  impulse  along  the  chain,  there  may  be  energy  losses  and 
we  consider  two  cases  described  in  the  following  subsections. 


3.1.1  Lossless  STC  Hard-Sphere  Approximation 

Ignoring  any  energy  loss  during  a  collision,  we  perform  an  STC  hard-sphere  approximation  for  a 
TC.  Masses  and  radii  are  expressed  as 

rM=ri-riq  =  (\-q)ri=sri, 


m.  =  pv.  = 


4 

I 


7irf  p  =q  rf. 


(1) 


and 


=  rjrf+l  =  qi  rf , 


(2) 


where  s=  1  -  q.  Evaluating  the  conservation  of  momentum  with  a  single  prime  denoting  post 
collision  values  and  the  initial  condition  that  the  /  + 1  particle  is  stationary  before  a  collision 
( v;+1  =0),  all  77  cancel  and  we  obtain 


mivi  +  mMvM  =  mlvl  +  mMvM 

3  3  '  3  3  f 

rivi  =  rivi  +sr:vM  , 


'  3 

V,  =  V;  +6*  V 


+1 


(3) 
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where  equations  1  and  2  have  been  used.  Following  the  same  procedure  for  the  conservation  of 
energy  while  ignoring  the  factor  of  one -half  yields 


v?=v;2+^. 


(4) 


Letting  A  =  £3v'M,  we  can  rewrite  equation  3  in  terms  of  v\  and  substitute  the  resulting 
expression  into  equation  4, 

v?=(vf—  A)2 +  Av'm 
=  v2  -2Avi  +  A2  +  AVm  , 

2v.  =  A  +  v'm 

y;+1  2 

V,.  1  +  £ 3 


Note  that  for  one  collision, 
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For  N  particles  there  will  be  N- 1  collisions,  each  of  which  has  the  ratio  in  equation  6.  Therefore, 
the  normalized  KE,  KEN ,  for  the  lossless  STC  hard-sphere  approximation  is  given  as 


ken 


4(1-?)3 
1  +  (1-?)3 


(V) 


3.1.2  Lossy  STC  Hard-Sphere  Approximation 

The  same  approximation  can  be  performed  with  some  amount  of  energy  loss,  E , ,  included  such 
that  the  system  is  still  conservative  and  may  represent  a  better  approximation.  Consequently,  the 
momentum  equation  3  is  unchanged,  but  equation  4  becomes 


v2=v'2+s\'2x  +  EL. 


(8) 


One  then  obtains  a  more  complicated  expression  replacing  equation  5: 


vl+i  _  gVm 

v,  l  +  e 3 

We  can  make  the  substitution,  EL  ocv;v'+1or  EL  =  E ,  v;v' , , ,  where  EL  is  the  constant  of 
proportionality.  This  adjustment  yields 


(9) 
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v;+i_  2 e3-EL 
v,  ^3(l  +  ^3)  ’ 

The  corresponding  result  for  the  normalized  KE  for  N  particles  is 


KEn=\ 


2(1  -qf~EL 


(1  -qf  l  + 


N-\ 


(10) 


In  the  limit  EL=  0 ,  equation  10  reduces  to  the  lossless  case,  equation  7,  as  one  would  expect. 
Note  that  results  are  independent  of  initial  velocity  and  size  of  the  grains. 

3.1.3  KE  Parameter  Space  for  STC  Hard-Spheres 

Figure  2  highlights  the  behavior  of  equation  10  for  0  <  q  <  0. 1 , 3  <  N  <  20  and  selected  EL .  The 
tapering  q  resembles  a  sigmoid  or  half-gaussian  and  is  stretched  to  infinity  for  small  N.  In  the 
limits  that  q  =  0  or  N=  1  (not  shown),  KE  is  unity.  For  the  lossless,  monodisperse  chain  (panel 
(a),  q  =  0),  that  energy  is  completely  transferred  regardless  of  the  number  of  spheres. 

Interestingly,  if  the  initial  velocity  is  supplied  to  the  smaller  end  of  the  TC,  one  also  observes 
shock  absorption  similar  to  the  system  in  figure  1,  albeit  with  less  efficiency.  This  is 
accomplished  by  adjusting  the  definition  of  tapering.  In  this  case,  subsequent  particles  are 
growing  in  size,  i.e.,  ri+l  =  (1  +  q)r,.  Equations  7  and  10  are  modified  accordingly,  and  the  result 
for  a  lossless  system  is  illustrated  in  figure  3.  Apparently,  both  configurations  mitigate  a 
propagating  pulse. 

3.2  The  Decorated  Tapered  Chain 

The  second  and  last  TC  geometry  that  we  consider  is  the  decorated  tapered  chain  (DTC) 

(figure  4).  This  can  be  assembled  from  the  STC  by  introducing  a  single-sized  interstitial  grain  of 
radius  reven  between  every  member  of  an  STC.  We  constrain  the  system  to  an  odd  number  of 
particles  such  that  the  inserted  grains  are  not  at  the  boundary.  Additionally,  we  presume  that 
these  interstitial  grains  will  be  equal  to  or  smaller  than  the  smallest  member  ( rN )  of  the  tapered 
part  of  the  chain;  therefore,  reven  =  frN,  where  0  <  /  <1.0 — although  the  flexibility  is  already 
built  in  for /  to  be  any  size. 
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Figure  2.  KE  ( N,q,EL )  parameter  space  for  the  STC  hard-sphere  approximation.  N  varies  from  3  to  25  and  q  from 
0%  to  20%. 
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q(%) 


Figure  3.  KE{N,q)  parameter  space  for  the  STC  hard-sphere  approximation  where  the  initial 
velocity  is  supplied  to  the  smaller  end  of  the  chain. 


Figure  4.  The  decorated  tapered  chain:  N=  13,  q  =  8%,/=  0.7,  L  =  80.7  mm,  and  r,  =  5  mm. 

It  is  immediately  clear  that  size  mismatch  between  neighboring  grains  is  a  function  of  position 
along  the  DTC.  This  is  in  stark  contrast  to  the  STC,  where  successive  grains  are  always  smaller 
(or  larger)  by  the  same  amount.  It  is  possible  then  to  have  DTC  chains  that  appear  to  resemble 
monodisperse  chains  for  only  a  portion  of  the  chain. 
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Again,  our  primary  interest  is  in  deriving  an  expression  for  the  normalized  KE, 
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(11) 


which  can  again  be  obtained  by  the  conservation  of  momentum  and  energy  where  collisions  are 
treated  as  independent.  Since  these  are  not  trivial  iterative  solutions  as  in  the  STC  hard-sphere 
approximation,  we  must  solve  for  successive  collisions  and  determine  the  overall  pattern.  We 
will  eventually  look  for  fonns  of  v'M/vi  and  then  generalize  for  N particles  or  A  - 1  collisions. 
First,  the  relationship  among  masses  and  radii  must  be  evaluated.  Given  A(odd)  particles  in  a 
DTC,  every  other  grain  will  reduce  in  size  by  q%  such  that  r.+2  =  srt where  e  =  (1  -  q) . 
Assembling  the  radii,  we  have 

r,  (12) 

rM  =frN  (!3) 


ri+ 1  =  n  -  qr  =  (1  -  q)r  =  sr.  (14) 

rM  =frN  (!5) 

ri+ 4  =  ri+ 2  -  qri+ 2  =  (!  -  q)ri+ 2  =Kri  (16) 

ri+5=frN  (17) 

ri+ 6  =  ri+ 4  -  qrM  =  (1  -  q)ri+4  =  ( 1 8) 


rN-i=frN  (!9) 

rN  =  fv-2  -  Vn-2  =  (1  -  4)rN-  2  =  sN-]>!2r  (20) 

Equations  12-14  imply  N  =  3,  equations  12-16  imply  N  =  5,  etc.  The  major  equations  for  radii 
are  therefore 


r  _  „(^-l)/2 
rN  ~  b  'i 


q, •+!),(, •+3),-,(JV-l) 


=  fs 


(21) 
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Recall  for  masses  that  mi  =  pVl  =  —  7tr?  p  =  rjrf . 


Note  that  since  r)  is  just  a  constant  and  will 


cancel  once  the  conservation  equations  are  put  into  use,  we  will  ignore  it  from  now  on.  This 
expression  for  mi  combined  with  equations  for  rt  provide  the  relations 
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(22) 


mN-X  *  VN- 1  =  Ari 
3(N-\)/2  3 

mN  oc  s  '  rt 


where  A  =  /3<s'3<v  l->/2.  We  may  now  use  these  to  evaluate  the  conservation  equations.  Beginning 
with  momentum  and  assuming  that  each  subsequent  particle  in  the  chain  begins  at  rest,  we  solve 
for  the  first  five  collisions.  Primes  and  double-primes  indicate  post-collision  states.  A  primed 
quantity  denotes  the  first  post-collision  state  of  a  sphere  which  serves  as  input  to  the  next 
collision.  To  keep  track  of  its  velocity  after  the  second  collision,  it  is  denoted  by  a  double-prime 
and  will  eventually  be  eliminated.  With  s  =  (1  -  q) ,  we  obtain 


mivi=miv'i+mMv'M 

-> 

vt  =  v'  +  Av'm 

(23) 

™i+yM  = 

:  mi+ivi+i  +  mi+2v'i+2 

->■ 

Av'M=Av”M+£yi+2 

(24) 

m,+2V'i+2 

=  mi+2Vi+2+mi+3V'i+3 

->■ 

£3<2=£3vl2  +  AVi+3 

(25) 

mi+3V'i+3 

=  mi+3v';+3  +  m[+4v'i+4 

-> 

Av'+3  =  Avy+s6v;+4 

(26) 

mi+yM 

= ml+yy + mi+y+5 

£6v'+4  =£6v:4+Av'+5 

(27) 

From  the  pattern  in  the  previous  fonnulas,  equation  23  can  be  rewritten  as  s\\  =  £°vj  +  Av'j+l. 
Evaluating  the  energy  conservation  equations  yields  the  same  form  as  equations  23-27  except 
velocities  are  squared: 


vf  =  v'2  +  Av 


i+i  ’ 


Av'2  -  Av"z  i  c. V2 

AVi+ 1  —  AVi+ 1  +  6  Vi+2  ’ 


3  1 2  3  ttl  .  a  f2 

£  Vi+2=£  VM+  Avm, 


(28) 

(29) 

(30) 


9 


(31) 

(32) 


A  tl  A  ft  1  .  fl 

Avi+3=Avi+3  +  £  v,.+4, 


.,6,  f2 


„6 


j  2 


S"VM  ~  £"V'i+4  +  Av'i+5 . 


We  can  combine  equations  23-27  and  28-32  to  eliminate  the  double-primed  terms  and  form  the 
y'  y' 

velocity  ratios:  — —  ,  — — ,  etc.  Beginning  with  equation  23,  we  isolate  v'  and  square  it  to 
v,.  v'., 

obtain  v'2  =  v2  -2Aviv'+1  +  A1v'l1 .  Next,  substitute  this  into  equation  28  and  rearrange  to  obtain 
v' 

— L±L.  This  is  then  repeated  for  equation  pairs  24  and  29,  25  and  30,  etc.,  to  obtain  the  following 

v,. 

ratios: 
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(34) 

(35) 

(36) 


where  we  insert  a  tenn  of  £°  in  equation  33.  With  our  goal  being  relation  1 1,  we  combine 
equations  33-36: 
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(38) 
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(39) 


After  some  observation,  the  ratio  can  be  put  into  closed  fonn  to  obtain 
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Turning  to  the  mass  ratios,  it  appears  that  most  tenns  cancel. 
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leading  to  the  simple  expression 

mN  _  3(JV-l)/2 

m1 


(41) 

(42) 


(43) 


We  can  now  identify  the  nonnalized  KE,  equation  1 1,  by  combining  equations  40  and  43  to  fonn 


KEn=(4As3/2) 


(JV-l) 


(7V-l)/2 
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with  A  =  /V(A  1)/2  and  £■=  (1  -£/).  This  is  plotted  in  figure  5,  where  it  is  clear  that  as/ 
decreases — so  that  inertial  mismatches  increase — the  normalized  KE  decays  rapidly.  Note  that 
the  colorscale  is  calibrated  for  each  subplot. 


4.  Numerical  Solution  to  the  Equations  of  Motion 


A  more  accurate  method  of  evaluating  the  dynamics  of  TCs  is  to  fonnulate  and  solve  the 
differential  equations  of  motion  (EOM).  Given  that  N  is  typically  small,  F  -  mizi  can  be 

solved,  allowing  us  to  easily  keep  track  of  velocities,  forces,  energies,  and  positions.  The 
problem  will  consist  of  TCs  barely  touching  in  their  initial  configuration  and  affixed  between 
rigid  walls  or  spheres  of  infinite  radius. 
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Figure  5.  KE  as  a  function  of  N,q  and  /  for  DTC  hard-spheres.  Note  that  the  colorscale  is  calibrated  for  each 
subplot. 


The  always-repulsive  potential  between  adjacent  grains  i  and  i+ 1  is  due  to  geometric  effects 
and  derived  in  Landau  and  Lifshitz  (23).  It  can  be  written  as 


5D\jR;+R 


RjRi+\  g 5/2 

uU+ 1  ui,i+l  i,i+l  - 


S5/ 2 


i+l 


(45) 


where  R. :  j  =  {/,  i  + 1,  •  •  • ,  N]  are  the  radii  and 

1  -  erf  1  -  cr 


D  =  - 
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(46) 
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Here,  8t  i+1  =  R:  +  RM  -  (zM  -  z;)  represents  the  overlap  between  successive  grains,  where  z .  is 
their  position.  Additionally,  aiM  has  been  defined  for  material  properties:  Ej ,  the  Young’s 
modulus  and  crj ,  the  Poisson’s  ratio;  and  radii,  i? . .  Note  also  that  j  can  refer  to  either  particle  i 
or  i  +  1 .  If  8i  i+x  <  0 ,  then  V  =  0  since  adjacent  grains  i  and  i  + 1  begin  to  lose  contact.  For  our 
particular  study,  all  materials  are  the  same  in  any  single  TC,  so  it  reduces  to 


D  = 


1  -  a 


2V 


(47) 


Under  the  Hertzian  potential,  equation  45,  a  squeezed  sphere  initially  appears  soft.  As 
compression  increases,  its  stillness  increases  substantially.  In  this  sense,  a  Hertz  potential  can  be 
thought  of  as  a  nonlinear  spring. 

At  the  boundaries,  the  enclosing  walls  are  represented  by  spheres  of  infinite  radius:  R0  and  RN+l. 
The  corresponding  potentials  then  for  grains  against  the  boundary  scale  as 

^(<L//,i)ccy/^  and  V{dNwall) kJr^. 


The  force  on  grain  i,  therefore,  is  the  sum  of  influences  from  its  neighboring  grains: 


mz 


dV 

dS 


(48) 


These  equations  are  solved  numerically  using  the  Velocity-Verlet  algorithm  (27).  The  original 
source  code  (5)  which  applies  to  the  STC  is  listed  in  appendix  A.  Modifications  to  this  for  the 
DTC  are  listed  in  appendix  B. 

Results  have  been  obtained  for  a  large  selection  of  chains  consisting  entirely  of  ThAfiV  or  SiC 
spheres.  Arbitrarily,  we  have  chosen  to  use  Ti6Al4V  when  restitutive  losses  are  ignored  ( co  =  0 ) 
and  SiC  otherwise.  The  following  material  properties  were  assumed  (reference  D  in 
equation  47). 


Tablet.  Material  properties  (28). 


Material 

P 

(mg/mm3) 

D 

(mm2/kN) 

Occurrence 

SiC 

3.2 

0.003266 

co  ^  0 

Ti6Al4V 

4.42 

0.01206 

co  =  0 

Note  that  the  fundamental  unit  of  force  in  our  simulations  is  the  kilonewton  and  an  initial 
velocity  of  0.01  mm/ps  is  applied  to  the  largest  end  of  each  chain. 

Simulations  were  performed  for  3  <  N  <  20  and  0  <  q,co  <0.1  over  a  system  time  of  1  ms  where 

o 

the  timestep  was  set  to  10  ps,  corresponding  to  10'  steps  in  the  integration  loop. 
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Restitution,  the  mechanism  for  introducing  energy  losses,  is  tuned  by  an  asymmetry  between  the 
(contact)  forces  of  loading  and  unloading  during  a  collision: 


F 


unloading 


F, 


=  1  -co 


loading 


(49) 


Therefore,  the  unloading  or  expansive  force  in  a  collision  is  some  fraction  of  the  loading  force 
and  perfectly  elastic  collisions  correspond  to  co  =  0. 

In  order  to  handle  the  hundreds  of  TC  calculations,  a  Practical  Extraction  and  Report  Language 
(PERL)  script  (appendix  C)  was  created  to  automate  the  process.  It  consists  of  nested  loops  that 
iterate  through  N  and  q  inserting  their  new  values  into  a  template  containing  the  C  code  which 
solves  the  EOM  of  the  TCs.  It  then  copies  this  to  a  directory  labeled  by  some  appropriate 
mnemonic.  This  was  done  so  that  when  visualizing  the  data  in  MATLAB,  it  would  be  easy  to 
change  directories  as  part  of  a  loop.  An  example  directory  structure  is, 

/w0/N2  0/N2  0ql0/  taperchainl .  cpp.  Data  files  are  created  for  every  grain  in  each  TC 
and  their  number  is  encapsulated  in  the  filename.  Appendix  D  contains  the  MATLAB  code  to 
generate  KE  surfaces. 

Recall  that  we  are  interested  in  forming  the  ratio  KEN  =KEoujKEin  ,  the  normalized  KE.  Here, 
KEout  is  the  first*  kinetic  energy  peak  felt  by  the  last  (smallest)  grain  and  KE in  is  a  constant 
determined  from  the  initial  velocity  and  density  of  the  first  grain.  Conversely,  Fout  is  the  first 
minimum  felt  by  the  last  grain  since  the  direction  of  a  negative  force  is  into  the  wall  or  force 
sensor.  The  algorithm  to  pick  out  the  first  turning  point  for  each  of  these  is  straightforward. 

Since  KE(t )  is  a  column  vector,  one  can  iterate  through  each  element  until  the  first  occurrence  of 
(i  + 1)<  i  is  true.  In  that  case,  i  represents  the  extremum.  For  force,  we  compare  against 
(/  + 1)  >  i .  The  only  data  files  then  that  are  of  interest  for  KE  ,.  are  those  that  represent  the 
dynamics  of  the  first  (KE  in)  and  last  particle  ( KE  out )  in  each  chain. 


4.1  STC 

Figure  6  plots  a  small  cross  section  of  the  results  for  several  chains  with  A  =  20,  co  =  0.05, 
andg  =  {0,  0.02,  0.04,  0.06,  0.08,  0.1}.  These  plots  show  the  kinetic  energy  spectrum  of  the 

tail  particle  as  a  function  of  time  for  the  first  half  of  the  simulation  (500  ps).  From  equation  49, 
each  collision  only  transmits  95%  of  the  force.  Note  the  decreasing  scale  of  KE  in  each 
underlying  panel.  One  might  expect  KE  to  increase  since  particle  velocity  is  increasing. 
However,  KE  scales  as  v2  and  m  oc  r3  oc  (1  -qf ,  so  the  latter  term  dominates,  and  consequently 
KE  decreases. 


In  the  top  panel  of  figure  6,  we  see  at  about  95  ps  that  the  tail  particle  in  a  monodisperse  chain 
first  receives  the  energy  bundle  and  undergoes  translational  motion.  Its  contact  force  with  the 


Normalization  is  discussed  further  in  appendix  E. 
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Figure  6.  KE  as  a  function  of  time  for  the  smallest  particle  in  a  chain  with  N  =  20,  CO  =  0.05,  and 

q  =  [0  :  0.02  :  0. 1].  Each  plot  is  evaluated  for  a  different  tapering  with  initial  KE  as  0.0838  J. 
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boundary  increases  as  it  slows  down  eventually  coming  to  rest  with  all  energy  stored  as  elastic 
potential  at  about  100  ps.  This  is  then  converted  back  to  kinetic  as  the  tail  particles  begin  to 
rebound. 

By  selecting  the  first  peak  in  each  of  these  subplots  as  well  as  the  many  other  chains  not  shown 
here,  we  can  form  the  nonnalized  KE  surfaces  and  observe  how  well  STCs  perform  as  energy 
absorbers. 

Pulses  in  monodisperse  chains  ( A>  15,  q  =  0)  represent  a  special  class  of  waves  known  as  a 
solitary  waves  which  are  remarkably  stable  and  have  been  studied  extensively  (2-7,  9,  11-12). 
For  TCs  (q  0),  the  solitary  wave  cannot  quite  fully  fonn  because  subsequent  grains  within  the 
wave  envelope  move  at  higher,  disparate  velocities  with  increasing  g — leading  to  greater 
dispersion  and  destructive  wave  broadening.  This  translational  symmetry  breaking  was  also 
reported  by  and  is  in  agreement  with  Nakagawa  et  al.  ( 4 ). 

4,1,1  STC  KE  and  Force  Parameter  Spaces:  KE  (N,q,co) ,  F  (N,q,co) 

Figure  7  highlights  the  numerical  results  for  KE(N,  q,  constant  co).  Note  that  it  would  be  very 
difficult  and  costly  to  produce  a  similar  plot  empirically.  Each  node  on  the  surface  represents  a 
different  TC  and  experiment.  There  appears  to  be  a  sigmoidal  and  exponential  dependence  on  q 
and  A,  respectively,  in  agreement  with  the  hard-sphere  approximation.  These  KE  surfaces 
represent  STC  chains  that  thermalize  more  than  half  the  incident  energy  introduced  into  the 
system.  For  example,  the  least  effective  shock  mitigating  geometry  that  we’ve  simulated  here — 
KE(o)  =  0,  N  =  3,  q  =  0) — reduces  the  output  KE  by  about  60%.  That  amount,  after  restitution  is 
accounted  for  in  cases  where  co  ^  0,  is  distributed  among  the  other  grains  as  kinetic  and  potential 
energy.  Note  again  that  results  are  independent  of  initial  grain  size.  For  monodisperse  chains 
(q  =  0)  in  figure  7a,  there  is  asymptotic  behavior  as  N  increases  which  is  a  result  of  energy 
propagating  as  solitary  waves,  when  N  is  of  sufficient  number  (about  15),  energy  can  propagate 
with  negligible  loss.  Additionally,  there  is  a  rapid  drop  in  KE(3  <  N<  5,  q  =  0,  co  =  0)  which  is 
likely  due  to  the  width  of  the  incident  pulse  being  longer  than  the  TC.  This  implies  that  wave 
reflection  begins  before  the  impulse  is  fully  applied. 

A  similar  surface  can  be  plotted  representing  the  normalized  force  for  various  TCs.  This  is 
visible  in  figure  8  and  similar  to  the  KE  plots.  Again,  the  functionality  of  N  here  is  essentially  a 
one  to  two  phase  exponential,  but  the  sigmoidal  nature  of  q  is  less  obvious  and  can  be 
approximated  linearly  throughout  the  parameter  space.  In  general,  the  magnitude  of  F  appears 
to  be  twice  that  of  KE . 

The  numerical  results  can  be  compared  to  the  hard-sphere  approximation  by  scaling  the  latter  to 
the  largest  numerical  solution — KE  (N  =  3,  q  =  0) — and  plotting  their  difference  (figure  9). 
Energy  loss  is  not  accounted  for  in  this  comparison.  The  distribution  is  due  to  differences  in  the 
interaction  potential  and  width  of  the  energy  bundle.  Therefore,  it  appears  that  for  smaller 
particles  (i.e.,  larger  q)  compressive  effects  and  width  of  the  energy  packet  becomes  increasingly 
important.  This  is  compounded  further  by  increasing  the  number  of  particles,  N . 
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Figure  7.  Numerical  solution  of  KE(N,q )  parameter  space  for  constant  CO.  Note  that  it  would  be  very  difficult 
and  costly  to  produce  a  similar  plot  empirically.  Each  node  on  the  surface  represents  a  different  TC  and 
experiment. 


17 


Figure  8.  Numerical  solution  of  F(N,q )  parameter  space  for  constant  (0. 
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Figure  9.  Difference  plot  of  the  lossless  KE  parameter  spaces  for  the  hard-sphere 

approximation  (figure  2a)  and  numerical  solution  (7a).  Note  that  the  azimuthal  view 
has  been  rotated  180°  for  clarity. 

4.1.2  STC  Mathematical  Model 

It  would  be  convenient  to  have  a  predictive  capability  for  at  least  a  selection  of  individual  TCs  as 
a  function  of  N,  q,  and  to.  We  therefore  propose  a  mathematical  model  for  the 
KE  =  KE(a>,q,N  =  20)  parameter  space  of  the  form 

KE{o),q)  =  AeBqC eDo}E ,  (50) 

which  corresponds  to  a  two-dimensional  Weibull  distribution.  One  can  obtain  the  sigmoidal  and 
exponential  behavior  of  q  and  co,  respectively,  when  the  exponents  are  restricted  to  C  >  1 
and  is  <1.  For  simplicity,  we  set  E  =  1  and  C  =  3/2  since  these  values  provide  the  general 

functionality.  The  coefficients  B  and  D  were  evaluated  using  4th-order  polynomial  fits  and  the 
scaling  coefficient  A  is  essentially  the  point  KE(N  =  20, co  =  q  =  0)  determined  in  the  simulation. 

It  turns  out  that  a  2nd-order  fit  was  not  sufficiently  robust  and  a  higher-order  fit  yielded  marginal 
gains  at  the  cost  of  mathematical  encumbrance. 

This  model  currently  lacks  the  rigor  sufficient  for  planar  behavior  in  the  limit  of  small  N .  It  is 
unclear  at  this  point  if  N  can  be  completely  decoupled  from  co  and  q  and  written  as  an  additional 
exponential  term.  In  all  likelihood,  the  coefficients  B  and  D  would  be  written  as  functions  of  N . 
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Ultimately,  it  would  be  useful  to  expand  equation  50  and  have  a  model  that  completely  describes 
the  normalized  KE  as  a  function  of  N,  q,  co,  and  A  where  the  latter  is  a  constant  external  loading 
or  precompression  of  the  chain. 

In  evaluating  the  fits,  we  find  that  KE(co,  q,  N  =  20)  is  described  quite  well  by 


KE(co,q)  =  0.35544 -e1 


(  [-1.5055-10~V+4.016-10-V-3.981-10~V+0.0147?-0.05435]?3/2) 
(  [4.144-10~V+1.9554(rV-0.03962?2 +0.028874-8.341]®) 


(51) 


It  is  compelling  to  rewrite  equation  5 1  in  the  more  suggestive  form 


KE(a>,q)  oc 


„9/2 


„7/2 


n^1  ef!1 


(52) 


(53) 


where  the  powers  of  q  have  been  adjusted  to  be  written  as  a  series  in  n,  the  interaction 
potential — equivalent  to  2.5  for  Hertzian  spheres.  The  occurrence  of  an  apparent  series  in  n  is 
quite  striking.  Missing  from  this  of  course  is  a  (n  -2)  or  q  2  tenn  which  requires  a  different  fit 
for  B  in  equation  50.  Note  also  the  mixed  term  qco  in  equation  5 1 .  This  is  indicative  of  a  many- 
body  effect;  one  cannot  have  restitution  in  a  single  particle  chain  since,  by  definition,  it  requires 
at  least  two  to  evaluate  the  loading  and  unloading. 

Figure  10  compares  the  numerical  results  and  the  model  for  the  case  of  N  =  20.  The  results  are 
practically  indistinguishable  except  for  minor  timestep  errors  in  the  simulation  which  appear  as 
noise  for  large  q.  At  these  values  for  N  and  q,  the  tail  particle  has  a  large  velocity  and  for 
timesteps  not  small  enough,  the  KE  plots  lose  their  smoothness.  The  algorithm  therefore  selects 
a  peak  close  to  what  the  real  extremum  would  be  in  the  limit  that  the  timestep  approaches  zero. 
As  such,  the  model  is  more  accurate. 


4.2  DTC 

In  this  section,  we  address  the  more  capable  DTC.  KE  time  spectra  is  excluded  since  an  analysis 
of  the  complicated  motion  is  outside  the  scope  of  this  report.  Instead,  let  us  focus  immediately 
on  the  KE  parameter  space  and  compare  with  the  STC. 

Figure  1 1  illustrates  the  intriguing  nature  of  the  DTC  KE(N,q,f )  surfaces  where  the 
dissimilarity  from  the  DTC  hard-sphere  approximation  is  quite  clear.  Immediately  obvious  is  a 
ripple  across  the  surface  for  large/ which  translates  towards  the  origin  as/  decreases.  The  effect 
vanishes  at  about/=  0.6  and  the  KE  surface  resembles  that  of  a  much-improved  STC  for 
/  <0.5.  We  propose  qualitatively  that  this  is  a  consequence  of  the  inertial  mismatch  changing 
as  a  function  of  position  in  the  DTC — in  contrast  to  the  STC  where  it  is  constant  along  the  chain 
and  the  wave  effect  is  not  observed. 
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Figure  10.  The  simulated  and  modeled  STC. 

In  the  side  panel  in  figure  1 1,  several  DTCs  are  shown.  In  the  top  case,/  =  1,  so  that  as  the 
incident  pulse  propagates  to  the  left,  the  chain  becomes  increasingly  monodisperse.  That  degree 
of  monodispersity  is  a  function  of  N  and  q  since  wave  amplitude  and  velocity  is  sensitive  to 
them.  If  N  is  large  enough,  then  whatever  remaining  energy  is  left  could  propagate  nearly 
lossless  like  a  solitary  wave.  As/ gets  smaller  (next  panel  down),  any  apparent  monodispersity 
begins  to  disappear.  As  one  moves  to  the  right  in  the  chain,  the  interstitial  grains  become  small 
compared  to  the  neighboring  grains.  As /  becomes  smaller  still,  the  interstitial  grain  is  small 
compared  to  all  grains  such  that  reducing/  has  little  effect — little  enough  so  that  the  wave 
behavior  on  the  KE  surface  is  not  observed.  Thus  there  appears  to  be  a  transition  so  that 
when  /  <0.5,  the  chain  will  never  appear  monodisperse  and  inertial  mismatches  can  always  be 

considered  large. 
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Figure  11.  DTC  numerical  solution  of  KE  (N,  q,f^)  and  sample  DTCs  for  various  f  (inset).  The  scale  is  the  same 
as  that  in  figure  7. 


5.  Experimental  Efforts 


Up  to  this  point,  most  experimental  work  reported  in  the  open  literature  has  focused  on 
monodisperse  chains  by  themselves  or  as  a  precursor  to  TCs.  Here  we  only  discuss  the  setup 
used  in  such  studies  as  well  as  the  prototype  for  our  intended  investigation  of  the  DTC. 

5.1  STC 

There  are  two  common  ways  of  performing  STC  experiments.  The  first  (11)  is  based  on  a 
Newton's  cradle  but  more  epic  in  extent  (figure  12).  In  this  particular  case,  there  is  a 
monodisperse  sequence  used  to  set  up  a  well-defined  solitary  wave  which  propagates  to  the  left. 
A  small  striking  sphere  on  the  far  right  is  applied  to  generate  the  initial  pulse. 
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Figure  12.  STC  experimental  setup  by  Job  et  al.  (72). 

The  second  common  STC  configuration  ( 4 )  consists  of  a  V-shaped  track  which  can  be  slanted  to 
ensure  proper  (centerline)  contact  of  the  tapered  spheres  (figure  13).  Here,  the  left  column 
highlights  a  monodisperse  chain  and  the  right  column  shows  the  same  experiment  but  for  a  TC — 
subsequent  rows  represent  later  time.  In  the  middle  row,  the  striker  (gold  sphere)  has  hit  the 
chain;  in  the  bottom  row,  the  tail  particle  in  the  monodisperse  chain  has  perforated  the  paper 
target  (black  hole  in  upper  left  of  picture),  while  for  the  TC  on  the  right,  it  merely  dents  and 
bounces  off  the  paper.  This  is  a  remarkable  demonstration  of  the  impulse  decimation  capability 
inherent  to  TCs. 

5.2  DTC 

Figures  14-16  highlight  the  constituents  and  current  housing  of  a  DTC  configuration. 

Within  the  assembly,  interstitial  grains  are  radially  constrained  along  the  axis  of  the  DTC  by  thin 
discs  with  small  holes  bored  in  them  (figure  16).  The  apparatus  can  be  tested  in  any  orientation; 
however,  when  vertical,  asymmetric  loading  due  to  gravity  will  become  a  factor.  Uniform 
precompression  enhances  the  shock  absorption  capability  of  TCs  and  will  be  discussed  in  the 
next  section.  In  figure  15,  we  see  that  precompression  of  the  DTC  prototype  can  be  controlled 
by  the  large  screw  at  the  top  of  the  assembly. 


23 


Source:  Picture  courtesy  of  Prof.  M.  Nakagawa  of  the  Colorado  School  of  Mines  and  Dr.  Juan  Agui  of  the  National  Aeronautics 
and  Space  Administration-Glenn  Research  Center. 

Figure  13.  STC  experimental  setup.  The  left  half  of  the  figure  represents  a  monodisperse  chain  while  the  right  a 
TC.  Lower  plots  denote  later  time. 
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Figure  14.  Preliminary  DTC  housing  chamber. 
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Figure  15.  DTC  constituents. 
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6.  Precompressed  Chains 


Apparently,  when  TCs  are  precompressed,  shock  absorption  increases  substantially  (J).  We 
measure  precompression,  A,  as  a  percentage  of  particle  overlap  in  the  initial  state.  The  EOM  in 
equation  48  are  modified  accordingly  and  become 


m,.z. 


^  K,„  K,  +*,  -  A-(z,  2 -aUtl  [«,+«,. 


Figure  17  highlights  the  impressive  results  for  KE(N  =  20,  A  =  0.03%)  as  tapering  and 
restitution  are  varied.  Precompression  in  a  DTC  could  therefore  provide  the  best  shock 
absorption  geometry  that  we  have  currently  considered. 


(54) 
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Figure  17.  KE  as  a  function  of  tapering  and  restitution  for 
constant  N  =  20  and  A  =  0.03%  . 


7.  TC  Armor  Panels 


We  have  seen  that  Hertzian  TCs  act  as  shock  absorbers.  In  particular,  recall  that  at  the  end  of 
section  3. 1 .3  we  noted  that  TCs  operate  regardless  of  a  pulse's  direction  through  the  chain.  Thus, 
to  maximize  efficiency,  we  can  place  neighboring  STCs  in  alternating  directions. 

For  them  to  become  realizable  in  armor  applications,  however,  one  should  evaluate  them  based 
on  performance  of  the  whole  system.  The  specific  absorbed  energy  (SAE) — which  has  units  of 
J/g — is  one  such  metric  and  commonly  used  by  ballisticians.  It  is  quantified  by  equation  55, 
while  the  volumetric  absorbed  energy  (VAE)  is  described  by  equation  56  and  has  units  of  J/cm  . 


KE.  - KE 

I  in  o. 

M 


x  KE0U, 

KE ... 


KE. 


M 


|l  -  KE (N,  q,a>)\KEin 
M 


|l  -  KE  (A,  q,  co)\ KEjn 
V 


(55) 

(56) 


28 


Here,  KE(N,q,a> )  is  the  normalized  KE  that  has  been  evaluated  in  the  numerical  simulations. 
Figure  18  illustrates  what  a  TC  armor  prototype  might  resemble  along  with  an  incident  flyer 
plate. 


Figure  18.  STC  armor  plate  panel. 


In  equations  55  and  56,  KE  jn  represents  the  energy  of  the  incoming  flyer  plate.  For  convenience, 
we  assume  that  all  materials  are  Ti6Al4V  ( p=  4.42  mg/mm3),  that  uf=  100  m/s,  t  =  24.5  mm, 
and  that  l=w  =  230  mm.  Since  m  =  pV  =  plwtf,  the  KE  of  the  flyer  plate,  KEin ,  is  28.7  kJ. 
Also,  in  both  cases,  we  set  KE(N,q,co)  =  KE (20, 0.1, 0.06)  «  0.05,  which  is  obtained  from  the 
numerical  results.  Thus,  the  sizes  of  the  end  spheres  in  such  a  chain  are  d20  =  10  mm  and 
dx  =  1.35  mm. 

In  order  to  evaluate  M  in  equation  55,  several  assumptions  are  made  and  then  later  relaxed  to 
determine  a  range  for  M  and  better  approximate  what  type  of  performance  a  TC  armor  panel 
might  achieve.  Figure  18  displays  the  relevant  variables  and  configuration  for  a  preliminary 
analysis,  where  t  =  95.8  mm  is  chosen  based  on  the  chain  parameters  used  in  picking 
KE(N,q,cd). 

Each  STC  would  require  some  type  of  mechanical  support  to  maintain  the  overall  geometry. 

One  solution  is  to  house  them  in  hollow  cones  as  is  displayed  in  the  diagram.  Conical  springs,  or 
Belville  washers,  are  another  possibility.  We  can  absorb  this  added  material  into  dx  as  an 
effective  diameter  which  we  set  to  c/,  =1.5  mm.  The  additional  dimensions,  w,l,  are  arbitrarily 
assigned  to  be  a  function  of  the  number  of  TCs:  l=w  =  20  (d20  +  <7,)=  230mm,  so  that  each  side 
is  spanned  by  40  TCs.  Therefore,  the  total  number  of  TCs  in  the  armored  panel  is  1600.  Finally, 
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St v St 2  represent  thicknesses  for  enclosing  plates  so  that  the  TCs  aren't  exposed.  This  amount  is 
arbitrary  but  chosen  to  be  small,  St]  =  St2  =  2.5  mm,  so  that  the  total  thickness  of  the  panel  is 
100.8  mm,  or  about  4  in. 

Now,  M  can  be  written  as  the  sum,  M  =  NSTCmSTC  +  2 ma  +  /nsupport  +  mm,  where  the  first  term  is 

the  product  of  the  mass  of  one  STC  and  the  total  number  of  them,  the  second  is  the  contribution 
from  the  enclosing  plates,  the  third  is  due  to  the  conical  housing  of  STCs,  and  the  last  is  the  filler 
material  which  we  ignore  as  being  air.  Using  conservative  values,  we  let  ma  =  (5t)lwp  =  0.58 
kg,  Support  =  1  -0  kg  so  that  M  =  13.63  kg. 

Since  equation  55  scales  with  uf,  the  velocity  of  the  flyer  plate,  a  wide  range  of  SAE  can  be 
determined.  If  uf  =  100  m/s,  the  SAE=2.  Let  us  assume  then  that  6.63  <M<  20.63  and 
100  <  uf  <  300  m/s  and  plot  the  SAE  surface  which  is  visible  in  figure  19a. 

The  VAE  can  be  calculated  by  dividing  this  quantity  by  the  density  of  the  material — since  we  are 
assuming  the  whole  panel  is  Ti6Al4V  (figure  19b).  Another  measure,  however,  would  be  to 
divide  by  the  volume  of  the  entire  panel  since  that  space  is  already  claimed  and  cannot  be 
occupied  by  anything  else.  In  dividing  by  the  density,  we  are  only  counting  that  portion  of  the 
armor  panel  contained  by  T^AL/V;  it  doesn't  include  the  air  fills.  Note  that  the  SAE  and  VAE 
have  been  determined  primarily  analytically. 


Figure  19.  SAE  and  VAE  as  a  function  of  M ,  the  total  STC  armor  panel  mass,  and  U,,  the  incident  flyer  plate 
velocity. 
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8.  Concluding  Remarks 


This  report  has  focused  on  TCs  and  their  inherent  ability  to  decimate  propagaing  impulses.  Due 
to  a  lack  of  dependence  on  system  size,  it  appears  that  TCs  can  provide  meter  to  submicron  scale 
shock  absorption.  TCs  take  well-defined  pulses,  such  as  shocks,  and  turn  them  into  noise.  The 
signal  is  spread  out  over  time  and  space  and  the  process  is  known  as  thermalizing  an  impulse  (3). 

Essentially,  the  shock  is  stretched  from  a  few  microseconds  to  several  hundred  with  its  extent 
increasing  from  several  grains  to  all  members  in  the  TC.  This  is  in  opposition  to  monodisperse 
chains  which  act  as  highly  efficient  shock  transmitters  due  to  the  resilient  nature  of  solitary 
waves.  In  fact,  Poschel  and  Brilliantov  have  calculated  (29)  the  optimal  transmission  of  KE  for 
TCs. 

It  turns  out  that  regardless  of  whether  an  impulse  travels  forwards  or  backwards  along  a  chain, 
both  directions  mitigate  a  propagating  pulse.  At  later  times,  this  likely  competes  with 
coordinated  motion  of  adjacent  grains  which  may  create  unpredictable  energy  or  force  spikes  for 
certain  configurations,  as  figure  7  (lower  panel)  demonstrates.  Antiparallel  orientations  of  TCs 
can  then  be  arranged  to  exploit  the  collective  effects  of  many  TCs  as  an  armor  panel.  The 
capabliity  of  such  a  panel  is  analytically  derived  and  yields  values  for  the  SAE  on  the  same  order 
10  -10"  J/g  as  other  technologies.  However,  such  a  determination  is  still  rather  suggestive 
without  a  full  scale  prototype.  TCs  may  be  housed  in  hollow  cones  or  Belville  washers  (conical 
springs)  for  added  absorption.  Additionally,  there  is  evidence  (30)  that  nitinol  (a  nickel-titanium 
alloy)  may  be  a  better  material  to  choose  for  the  spheres  and  housing. 

Hard-sphere  approximations  were  analytically  obtained  for  two  chain  geometries,  the  STC  and 
DTC.  This  approximation  simplifies  the  problem  to  one  where  the  impulse  delivered  is  a  single 
grain  in  width  and  elastic  effects  are  ignored.  For  the  STC,  this  approximation  is  comparable  to 
the  numerical  results.  However,  for  the  DTC,  the  differences  are  quite  large  and  consequently 
elastic  effects  play  a  major  role  in  the  system.  This  is  likely  due  to  the  increased  velocity  of 
members  in  the  chain  and  therefore  the  number  of  collisions. 

The  authors  see  ample  opportunity  for  study  along  several  lines  to  move  this  technology  forward 
and  rate  it  as  a  useful  technology  in  blast  mitigation  system.  Quasistatic  compression  to 
plasticity  and  failure,  ALEGRA  simulations,  and  measuring  the  effects  of  many  TCs  working  in 
cooperation  are  areas  we  are  beginning  to  investigate. 
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Appendix  A.  Simple  Tapered  Chain  Code 1 


This  appendix  appears  in  its  original  form,  without  editorial  change. 

1  Pfannes,  J.  Energy  Propagation  in  Granular  Chains.  M.S.  Thesis,  State  University  of  New  York,  Buffalo,  NY,  May  2003. 
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/*  PROGRAM  taperchain . cpp 

This  program  consideres  an  one  dimensional  chain  of  spheres  that 
shrink  succesively  in  radius  ("tapered  chain") .  Initially  the  spheres 
are  barely  in  contact,  i.  e.  they  just  touch  each  other  and  are  not 
compressed  (zero  loading) . 

The  chain  ends  at  both  edges  at  fixed  walls. 

The  program  calculates  the  interaction  of  the  system  once  disturbed 
by  an  instantaneous  (delta)  impulse  exerted  on  one  end  of  the  chain. 
Restitution  both  between  the  spheres  and  between  the  edge  spheres 
and  the  corresponding  wall  can  be  introduced. 

The  program  does  not  consider  gravity. 

The  EOM  are  solved  with  the  Velocity  Verlet  algorithm. 

scale  of  problem:  mm-mg-musec  (mimimu) 
in  this  scale  unit  of  force:  1000  N 
in  this  scale  unit  of  energy:  1  J 

*/ 

#include  <cmath> 

#include  <iostream> 

#include  <fstream> 

#include  <cstdlib> 

#include  <string> 

#include  <sstream> 

using  namespace  std; 

/**************************  ALTERABLE  PARAMETER:  ****************************/ 
const  int  nptles=20;  //  total  number  of  particles 

const  double  rho=3.2  /*  SiC  (mg/mm^3)  */,  D=0 . 00326603139013  /*  (mm^2/N)  */; 

const  double  rlarge  =5.0;  //  (radius  of  large  ptle  (mm)) 

const  double  q  =  0.0;  //  (tapering  factor  (%) ) 

const  double  xn  =2.5;  //  (exponent  in  potential) 

const  double  dt  =  0.00001;  //  (timestepwidth  (musec) ) 

const  unsigned  int  nsteps  =  100000000;  //  (#  steps  integration  loop) 

const  int  idiagp  =  20000;  //  (stepwidth  diagnostics) 

const  int  idump  =  20000;  //  (stepwidth  dump) 

const  double  vlin  =  0.0;  //  (initial  v  small  ptle  (mm/musec) ) 

const  double  vnin  =  -0.01;  //  initial  v  large  ptle  (mm/musec)) 

const  double  epsilon  =1.0;  //  ( (1  -  restitution  factor)  all  ptles) 

/****************************************************************************/ 

ofstream  readme (" taperchain . readme ") ;  //  global  scope  fcts 

of stream  Energylmpulse (" taperchain . Enelmp" ) ; 

void  radii  (double  rlocal [] )  { 

rlocal [nptles-1]  =  rlarge; 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles-1;  i++) 
rlocal  [i]  =  rlarge; 

else 

for  (int  i  =  2;  i  <  nptles+1;  i++) 
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rlocal [nptles-i]  =  rlocal [nptles-i+1]  *  (1  -  q*0.01); 

} 

void  masses  (double  r[],  double  masslocal [] )  { 

const  double  pi  =  4  *  atan(l.O); 

const  double  masslarge  =  (4. 0/3.0)  *  pi  *  rlarge*rlarge*rlarge  *  rho; 
masslocal [nptles-1]  =  masslarge; 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles-1;  i++) 
masslocal  [i]  =  masslarge; 

else 

for  (int  i  =  0;  i  <  nptles-1;  i++) 

masslocal  [i]  =  r  [i] *r  [i] *r [i]  *  masslarge  /  (rlarge*rlarge*rlarge)  ; 

} 

void  strenghtfac  (double  r[],  double  alocal [] )  { 

alocal  [0]  =  (2.0  /  (5.0  *  D) )  *  (sqrt (r [0] ) ) ; 
alocal [nptles]  =  (2.0  /  (5.0  *  D) )  *  (sqrt (r [nptles-1] )) ; 
if  (q  =  =  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  1;  i  <  nptles ;  i++) 

alocal  [i]  =  (2.0  /  (5.0 *D) )  *  (sqrt (0 . 5*rlarge) )  ; 

else 

for  (int  i  =  1;  i  <  nptles;  i++) 

alocal [i]  =  (2.0  /  (5.0  *  D) )  *  (sqrt ( (r [i] *r [i-1] ) / (r [i] +r [i-1] ) ) ) ; 

} 

//  initialpos  prints  absolute  initial  positions,  not  for  calculations 
void  initialpos  (double  r[],  double  xlnitiallocal [] )  { 

if  (q  ==  0)  //  avoid  roundoff  errors  w/out  tapering 

for  (int  i  =  0;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  (2.0*  (i  +  1)  -  1)  *  rlarge; 

else  { 

xlnitiallocal [0]  =  r[0]; 

for  (int  i  =  1;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  xlnitiallocal [i-1]  +  r[i-l]  +  r[i]; 


//  absolutpos  prints  absolute  positions  to  ptle  files,  not  for  calculations 
void  absolutpos  (double  r[],  double  x[], 

double  xlnitial [] ,  double  xAbsolutlocal [] )  { 

for  (int  i  =  0;  i  <  nptles;  i++) 

xAbsolutlocal [i]  =  xlnitial [i]  +  x[i]; 

} 

void  computeAccelerations  (double  x[],  double  a[],  double  r[], 

double  acc[],  double  overbef ore [] , 
double  mass [] ,  doubles  pot)  { 

pot  =0.0;  //  every  call  calculates  new  pot  contributions 

for  (int  i  =  0;  i  <  nptles;  i++)  //  zeroing  all  acc  in  every  call 

acc[i]  =0.0;  //  (=  every  timestep) 

/*******  potential/f orce  between  neighboring  ptles  *****************/ 
for  (int  i  =  0;  i  <  nptles-1;  i++)  { 

if  (x[i]  >  x[i+l])  {  //  only  when  overlap 
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double  over  =  x[i]  -  x[i+l] ; 


double  overnml  =  pow(over,  (xn  - 
pot  +=  over  *  overnml  *  a[i  +  l]  ; 

1.0)  ) 

double  forceBetw  =  a[i+l]  *  xn  * 

overnml ; 

double  forceFactor; 

if  (overbef ore [i+1]  <  over) 

// 

when  compressing 

forceFactor  =  1.0; 

else  forceFactor  =  epsilon; 

// 

when  decompressing 

forceBetw  *=  forceFactor; 

// 

dim  acc 

:  force 

acc[i]  -=  forceBetw; 

// 

sign ( - ) 

:  towards  smaller 

acc  [i  +  1]  +=  forceBetw; 

// 

sign (+) 

:  towards  larger 

overbef ore [i+1]  =  over; 

} 

else  overbef ore [i+1]  =  0.0; 

// 

update 

for  next  timestep 

// 

reset  when  no  overlap 

/**  potential/force  between  fixed  wall  (small,  x=0)  <->  small  ptle  **/ 
if  (x  [0]  <  0)  { 

double  over  =  -  x[0]  ; 

double  overnml  =  pow(over,  (xn  -  1.0)); 

pot  +  =  over  *  overnml  *  a[0] ; 

double  forceSmall  =  a[0]  *  xn  *  overnml; 

double  forceFactor; 
if  (overbef ore [0]  <  over) 
forceFactor  =  1.0; 
else  forceFactor  =  epsilon; 

forceSmall  *=  forceFactor; 

acc[0]  +=  forceSmall; 

overbef ore [0]  =  over; 

} 

else  overbef ore  [0]  =  0.0; 

/***  potential/force  between  fixed  wall  (large)  <->  large  ptle  ******/ 
if  (x[nptles-l]  >  0)  { 

double  over  =  x[nptles-l]; 

double  overnml  =  pow(over,  (xn  -  1.0)); 

pot  +=  over  *  overnml  *  a[nptles] ; 

double  forceLarge  =  a [nptles]  *  xn  *  overnml; 

double  forceFactor; 
if  (overbefore [nptles]  <  over) 
forceFactor  =  1.0; 
else  forceFactor  =  epsilon; 

forceLarge  *=  forceFactor; 

acc [nptles-1]  -=  forceLarge; 

overbefore [nptles]  =  over; 
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} 

else  overbefore [nptles]  =  0.0; 

/*****  real  dim  of  acc :  division  by  mass  **********/ 
for  (int  i  =  0;  i  <  nptles ;  i++) 
acc  [i]  /=  mass [i]  ; 

} 

void  velocityVerletStep  (double  x[],  double  v[],  double  acc[], 

double  a[],  double  r[],  double  overbefore  []  , 
double  mass [] ,  doubles  pot)  { 

for  (int  j  =0;  j  <  nptles;  j++)  { 

x[j]  +=  v[j]  *  dt  +  0 . 5  *  acc[j]  *  dt*dt; 
v[j]  +=  0.5  *  acc[j]  *  dt; 

} 

computeAccelerations  (x,  a,  r,  acc,  overbefore,  mass,  pot) ; 

for  (int  j  =0;  j  <  nptles;  j++) 
v[j]  +=  0.5  *  acc[j]  *  dt; 

} 

void  ptleHeader  (of stream*  print,  int  k)  { 

(*  print)  <<  "#  ptle  "  <<  k+1  <<  time  (musec) "  <<  ' \t '  <<  "x  (mm)" 

<<  ' \t '  <<  "v  (mm/musec) "  <<  ' \t '  <<  "a  (mm/musecA2) " 

<<  ' \t ’  <<  "kin.  E.  (J) "  <<  ' \t ’  <<  "f  (kN) "  <<  ' \t ' 

<<  "impulse  (mg*mm/musec) "  <<  ' \t '  <<  "xRelative  (mm)"  <<  '\n'; 


void  dumpData  (double  t,  double  mass [] ,  double  v[],  double  acc[], 

double  r[],  double  x[],  //  scope  absolutpos 

double  xlnitial [] ,  double  xAbsolut [] , 
of stream*  print,  int  k)  { 

double  keDumpPt [nptles] ;  //  new  arrays  for  dumping  data 

double  vDumpPt [nptles] ;  //  since  arrays  pass  by  argument 

double  accDumpPt [nptles] ;  //  dump  data  manipulated 

double  xDumpPt [nptles] ; 

keDumpPt [k]  =  0.5  *  mass [k]  *  v [k] *v [k] ; 

if  (keDumpPt [k]  <  1.0e-20)  keDumpPt [k]  =  0.0;  //  set  small  values  to  zero 
vDumpPt [k]  =  v [k] ; 

if  (vDumpPt [k]  <  1.0e-20  &&  vDumpPt [k]  >  -1.0e-20)  vDumpPt [k]  =  0.0; 

accDumpPt [k]  =  acc [k] ; 

if  (accDumpPt [k]  <  1.0e-20  &&  accDumpPt [k]  >  -1.0e-20)  accDumpPt [k]  =  0.0; 

xDumpPt [k]  =  x [k] ; 

if  (xDumpPt [k]  <  1.0e-20  &&  xDumpPt [k]  >  -1.0e-20)  xDumpPt [k]  =  0.0; 

absolutpos  (r,  x,  xlnitial,  xAbsolut);  //  calculate  absolute  pos . 

(*  print)  <<  t  <<  ' \t '  <<  xAbsolut [k]  <<  ' \t '  <<  vDumpPt [k]  <<  ' \t 1 

<<  accDumpPt [k]  <<  ' \t 1  <<  keDumpPt [k]  <<  ' \t ' 

<<  accDumpPt [k] *mass [k]  <<  ' \t ' 
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<<  mass [k] *vDumpPt [k]  <<  ' \t '  <<  xDumpPt [k]  <<  1  \n'; 

} 

void  dumpEnergylmpulse  (double  t,  double  kelocal,  double  telocal,  double  pot, 

double  ptotallocal,  double  mass [] ,  double  v [ ] )  { 

kelocal  =  0.0; 
ptotallocal  =  0.0; 

double  absptotallocal  =0.0;  //  scope  only  within  function 

for  (int  j  =0;  j  <  nptles;  j++)  { 
kelocal  +=  masstj]  *  v[j]*v[j]; 
ptotallocal  +=  masstj]  *  v[j]; 
absptotallocal  +=  masstj]  *  abs(v[j]); 

} 

kelocal  *=  0.5; 
telocal  =  kelocal  +  pot; 
double  potDump  =  pot ; 

if  (kelocal  <  1.0e-20)  kelocal  =0.0;  //  set  very  small  values  to  zero 

if  (ptotallocal  <  1.0e-20  &&  ptotallocal  >  -1.0e-20)  ptotallocal  =  0.0; 
if  (telocal  <  1.0e-20)  telocal  =  0.0; 
if  (potDump  <  1.0e-20)  potDump  =  0.0; 
if  (absptotallocal  <  1.0e-20)  absptotallocal  =  0.0; 

Energylmpulse .precision (16)  ; 

Energylmpulse  <<  t  <<  ' \t '  <<  kelocal  <<  ' \t '  <<  potDump  <<  ' \t ' 

<<  telocal  <<  ' \t '  <<  absptotallocal  <<  ' \t ' 

<<  ptotallocal  <<  '\n'; 


void  readmelnfo  (double  ke,  double  pot,  double  te,  double  ptotal, 
double  kelin,  double  kenin,  double  plin,  double  pnin, 
double  r[],  double  xlnitial  []  ,  double  mass  []  , 
double  a  []  )  { 

readme  <<  1 \t 1  <<  ***  TAPERCHAIN  within  walls  ***  (-:" 

<<  ' \n '  <<  ' \n ' ; 


readme 

<< 

"parameter  of  this  run:  "  << 

1  \n ' 

<< 

1  \n ' 

/ 

readme 

<< 

"total  number  of  particles: 

"  << 

' \t 

'  << 

' \t  ’ 

<< 

nptles  << 

'  \n 

readme 

<< 

"density  of  particles  (mg/mm 

A3  )  : 

"  <<  ' \t 

1  << 

rho 

<<  1 \n' ; 

readme 

<< 

"quantity  D  of  particles  (mm 

^2/N) 

.  II 

<<  ' 

\t ' 

<<  D 

<<  1 \n' ; 

readme 

<< 

"radius  of  large  ptle  (mm) : 

"  << 

:\t 

'  << 

' \t  ’ 

<< 

rlarge  << 

'  \n 

readme 

<< 

"tapering  factor  (%) :  "  <<  ' 

\t '  <<  ' 

\t '  <<  ' \t '  <<  q  <<  ' \n 

1  . 

/ 

readme 

<< 

"exponent  in  potential:  "  << 

' \t  ’ 

<< 

' \t  ’ 

<< 

■  \t  - 

<  <  xn  <  < 

'  \n 

readme 

<< 

"timestepwidth  (musec) :  "  << 

' \t  ’ 

<< 

' \t  ’ 

<< 

■  \t  - 

<<  dt  << 

'  \n 

readme 

<< 

"#  steps  integration  loop:  " 

<<  1 

\t ' 

<<  1 

\t  ■ 

<<  nsteps  <<  1 

\n ' 

readme 

<< 

"stepwidth  diagnostics:  "  << 

' \t  ’ 

<< 

' \t  ’ 

<< 

' \t 

1  <<  "every  " 

<<  idiagp  <<  "  timesteps"  <<  '  \n'; 


readme 

<< 

"stepwidth 

dump : 

"  << 

' \t '  <<  ’ \t ' 

<< 

' \t  ’ 

<< 

"every  " 

<< 

idump  <<  " 

timesteps"  << 

1  \n '  ; 

readme 

<< 

"initial 

V 

small 

ptle 

(mm/musec) : 

"  << 

'  \f 

<< 

vlin  << 

1  \n 

readme 

<< 

"initial 

V 

large 

ptle 

(mm/musec) : 

"  << 

' \t ' 

<< 

vnin  << 

1  \n 

readme  <<  "restitution  factor  for  all  ptles:  "  <<  ' \t '  <<  1-epsilon  <<  '\n' 
<<  ' \n' ; 

readme  <<  "total  length  of  run  (musec) :  "  <<  ' \t '  <<  dt  *  nsteps  <<  '\n'; 

readme  <<  "total  rows  recorded  for  .Enelmp  file:  "  <<  ' \t '  <<  ' \t ' 
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<<  nsteps/idiagp  +1  <<  ' \n'; 

readme  <<  "total  rows  recorded  for  particle  files:  "  <<  ' \t ' 

<<  nsteps/idump  +1  <<  ' \n'  <<  1 \n'; 
readme  <<  "Initial  system  info  (t=0) :  "  <<  '\n'  <<  ' \n'; 
readme  <<  "kin.  E.  (J) "  <<  ' \t 1  <<  "pot.  E.  (J) "  <<  ' \t 1 

<<  "tot.  E.  (J) "  <<  ' \t '  <<  "total  impulse  (mg*mm/musec) "  <<  ' \n 
readme  <<  ke  <<  ' \t '  <<  pot  <<  ' \t '  <<  te  <<  ' \t '  <<  ptotal 
<<  ' \n '  <<  ' \n ' ; 


readme  <<  "kin.  E.  of  small  particle  (J) :  "  <<  ' \t ' 

<<  ' \n' ; 

<< 

' \t ' 

<<  kelin 

readme  <<  "kin.  E.  of  large  particle  (J) :  "  <<  ' \t ' 

<<  ' \n' ; 

<< 

1  \t ' 

<<  kenin 

readme  <<  "impulse  of  small  particle  (mg*mm/musec) :  " 

<<  pi in  <<  ' \n ' ; 

<< 

1  \t ' 

readme  <<  "impulse  of  large  particle  (mg*mm/musec) :  " 

<< 

' \t ' 

<<  pnin  <<  '\n'  <<  '  \n'; 
readme  <<  "particle  radii  (mm) :  "  <<  '  \n'; 

for  (int  i=0 ;  i  <  nptles;  i++)  { 
readme  <<  r[i]  <<  ' \t '  ; 

} 

readme  <<  '\n'  <<  '  \n'; 

readme  <<  "initial  particle  positions  (mm) :  "  <<  '  \n'; 
for  (int  i=0;  i  <  nptles ;  i++)  { 
readme  <<  xlnitial [i]  <<  ' \t '  ; 

} 

readme  <<  '\n'  <<  '  \n'; 

readme  <<  "total  length  of  one  dimensional  alignment  (mm) :  "  <<  ' \ t ' 

<<  xlnitial [nptles-1]  +  r[nptles-l]  <<  '\n'  <<  '  \n'; 
readme  <<  "particle  masses  (mg) :  "  <<  '  \n'; 
for  (int  i=0;  i  <  nptles;  i++)  { 
readme  <<  mass [i]  <<  ' \t '  ; 

} 

readme  <<  '\n'  <<  '  \n'; 

readme  <<  "particle  interaction  strenghts  (0 . 0316*N/mm^ (3/2) ) :  "  <<  ' \n 
for  (int  i=0 ;  i  <  nptles+1;  i++)  { 
readme  <<  a[i]  <<  '  \t '  ; 

} 

readme  <<  ' \n ' ; 

} 


int  main  (  )  { 

double  r [nptles],  x [nptles],  xAbsolut [nptles] ,  xlnitial [nptles] , 
v [nptles],  acc [nptles]  ,  mass [nptles]  ,  a [nptles  +  1], 
overbef ore [nptles+1] ,  pot ; 

/********************  functions  *****************************/ 

radii  (r) ; 

masses  (r,  mass) ; 

strenghtfac  (r,  a) ; 

initialpos  (r,  xlnitial) ; 

/**********************  output  files  ****************************/ 

/******  ptle-files  *********/ 

ofstream  print  [nptles] ; 

for  (int  k  =  0;  k  <  nptles;  k++)  { 
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string  filename; 
ostringstream  buffer; 

buffer  <<  " taperchain_"  <<  k+1  <<  ".dat"; 
filename  =  buf f er . str ( ) ; 

print [k] . open (filename . c  str());  //  convert  string  to  char 

} 

for  (int  k  =  0;  k  <  nptles;  k++)  //  header  for  particles 

ptleHeader (&  print [k] ,  k) ; 
j ***************************  j 

//  header  for  .Enelmp  file 

Energylmpulse  <<  "#  time"  <<  ' \t 1  <<  "kin.  E.  (J) "  <<  1 \t 1 


<< 

' \t ’  << 

"pot.  E.  (J) "  << 

' \t  ’ 

<< 

'  \t 

<< 

"total  E 

(J) "  <<  ' \t '  << 

•\t ' 

<< 

"  (total 

imp.) 1  (mg*mm/musec) " 

<<  1 

\t  ■ 

<< 

"total  imp.  (mg*mm/musec) " 

<< 

'  \n ' 

/ 

/*****************************************************************/ 

for  (int  i  =  0;  i  <  nptles ;  i++)  {  //  zeroing 

x[i]  =  0.0;  //  relative  particle  positions  for  calculation 

v [ i ]  =  0.0; 
acc[i]  =  0.0; 
overbef ore  [i]  =  0.0; 

} 

overbef ore [nptles]  =  0.0; 
double  t  =  0.0; 

v[0]  =  vlin;  //  mind  special  input  data 

v[nptles-l]  =  vnin; 

/****************  initial  energy  info  *********************************/ 
double  kelin  =  0.5  *  mass [0]  *  v[0]*v[0];  //  initial  kE  of  edge  ptles 

double  kenin  =  0.5  *  mass [nptles-1]  *  v [nptles-1] *v [nptles-1] ; 
double  plin  =  mass [0]  *  v[0];  //  initial  impulse  of  edge  particles 

double  pnin  =  mass [nptles-1]  *  v [nptles-1] ; 

computeAccelerations  (x,  a,  r,  acc,  overbefore,  mass,  pot) ; 

//  call  for  initial  potential  energy 

double  ke  =  0 . 0 ;  //  checking  for  initial  system  energy 

double  ptotal  =  0.0; 

for  (int  i  =  0;  i  <  nptles;  i++)  { 

acc[i]  =0.0;  //  reset  acc  for  verlet  in  case  of  initial  overlap 

ke  +=  mass  [i]  *  v[i]*v[i]; 
ptotal  +=  mass [i]  *  v[i]; 

} 

ke  *  =  0.5; 

double  te  =  pot  +  ke;  //  initial  total  system  energy 

readmelnfo  (ke,  pot,  te,  ptotal,  kelin,  kenin,  plin,  pnin, 
r,  xlnitial,  mass,  a) ; 

/*******************  begin  of  timestep  loop  *********************/ 
for  (unsigned  int  i  =  0;  i  <  nsteps+1;  i++)  { 
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t  =  i  *  dt  ; 

velocityVerletStep  (x,  v,  acc,  a,  r,  overbefore,  mass,  pot) ; 

/***********  check  system  energy,  plot  data  *******************/ 
if  ( (i  %  idiagp)  ==  0)  { 

dumpEnergylmpulse (t ,  ke,  te,  pot,  ptotal,  mass,  v)  ; 

} 

/*********  check  particle  energy,  plot  data  *******************/ 
if  ( (i  %  idump)  ==  0)  { 

for  (int  k  =  0;  k  <  nptles;  k++)  { 

dumpData(t,  mass,  v,  acc,  //  scope  dumpData 

r,  x,  xlnitial,  xAbsolut,  //  scope  absolutpos 

Sc  print  [k]  ,  k)  ; 


}/********************  encj  0f  timestep  loop  *********************/ 

readme . close ( ) ; 

Energylmpulse . close ( ) ; 

for  (int  k  =  0;  k  <  nptles ;  k++)  //  closing  particle  files 

print [k] . close ( ) ; 

} 
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Appendix  B.  Decorated  Tapered  Chain  Modifications 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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/**************************  ALTERABLE  PARAMETER:  ****************************/ 
int  nptles  =3;  //  total  number  of  particles 

const  double  f  =  0.9;  //  fractional  size  of  interstitial  grain  w.r.t  last 

grain 

const  double  rho=4.42;  //  TiAlV  (mg/mm^3) 

const  double  D  =  0.01206;  //  TiAlV  (mm^2/N) 


const 

double 

rlarge  =  5.0 

const 

double 

q  =  0 

-0; 

const 

double 

xn  = 

2.5; 

const 

double 

dt  = 

0 . 00001 

const 

unsigned  int 

nsteps 

const 

int  idiagp  = 

20000; 

const 

int  idump  = 

20000; 

const 

double 

vlin 

=  0.0; 

const 

double 

vnin 

=  -0.01 

const 

double 

epsilon  =  1. 

//  (radius  of  large  ptle  (mm)  ) 

//  (tapering  factor  (%)  ) 

//  (exponent  in  potential) 

//  (timestepwidth  (musec) ) 

100000000;  //  (#  steps  integration  loop) 

//  (stepwidth  diagnostics) 

//  (stepwidth  dump) 

//  (initial  v  small  ptle  (mm/musec) ) 

//  initial  v  large  ptle  (mm/musec) ) 

//  ( (1  -  restitution  factor)  all  ptles) 


/ 


ofstream  readme (" taperchain . readme ") ;  //  global  scope  fcts 

of stream  Energylmpulse (" taperchain . Enelmp" ) ; 


//  Generate  radii  and  masses  for  DTC 

void  spheres  (double  rlocal [] ,  double  masslocal [] )  { 


rlocal [nptles-1]  =  rlarge;  //  shift 
double  tapering  =  1  -  q*0.01; 
const  double  pi  =  4  *  atan(l.O); 
const  double  masslarge  =  (4. 0/3.0) 
masslocal [nptles-1]  =  masslarge; 
if  (q==0  &&  f==1.0)  // 

for  (int  i=0;  i<nptles-l;  i++) 
rlocal [i]  =  rlarge; 

masslocal  [i]  =  masslarge; 

else  if  (q==0)  { 

roundoff  errors  without  tapering 

for  (int  i=l;  i<nptles-l;  i=i+2) 
rlocal  [i]  =  rlarge*f; 
masslocal  [i]  =  (4. 0/3.0)  * 

} 

for  (int  i=0;  i<nptles-l;  i=i+2) 

grains 


;  everything  to  index  starting  at  zero 

pi  *  pow (rlarge, 3 )  *  rho; 

Monodisperse  in  DTC 

{ 

//  Quasi -Monodisperse :  Avoid 
{  //  Interstitial  grains 

pi  *  pow (rlocal [i] , 3 )  *  rho; 

{  //  Non-interstitial 


rlocal [i]  =  rlarge; 
masslocal  [i]  =  masslarge; 


II  non-monodisperse  chains 
else  { 

for  (int  i=l;  i<nptles-l;  i=i+2)  {  //  Find  radii  of 

interstitial  grains 

rlocal  [i]  =  f*pow (tapering,  (nptles-1) /2) *rlarge; 
masslocal  [i]  =  (4. 0/3.0)  *  pi  *  pow (rlocal  [i]  , 3 )  *  rho; 

for  (int  i=nptles-l;  i>=0;  i=i-2)  {  //Find  radii  of 

non- interstitital  grains 

rlocal  [i-2]  =  rlocal  [i]  *  tapering; 

masslocal  [i-2]  =  (4. 0/3.0)  *  pi  *  pow (rlocal [i-2]  , 3 )  *  rho; 


void  strenghtfac  (double  r[],  double  alocal [] )  { 

alocal  [0]  =  (2.0  /  (5.0  *  D) )  *  (sqrt (r [0]  ) )  ; 
alocal [nptles]  =  (2.0  /  (5.0  *  D) )  *  (sqrt (r [nptles-1])); 
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1.0) 


//  avoid  roundoff  errors  w/out 


if  (q  ==  0  ScSc  f  =  = 
tapering 

for  (int  i  =  1;  i  <  nptles;  i++) 

alocal  [i]  =  (2.0  /  (5.0*D))  *  (sqrt (0 . 5*rlarge) ) ; 

else 

for  (int  i  =  1;  i  <  nptles ;  i++) 

alocal  [i]  =  (2.0  /  (5.0  *  D) )  *  (sqrt ( (r [i] *r [i-1] ) / (r [i] +r [ i - 1] ) ) )  ; 

} 

//  initialpos  prints  absolute  initial  positions,  not  for  calculations 
void  initialpos  (double  r[],  double  xlnitiallocal [] )  { 

if  (q  ==  0  ScSc  f  ==  1.0  )  //  avoid  roundoff  errors  w/out 

tapering 

for  (int  i  =  0;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  (2.0*  (i  +  1)  -  1)  *  rlarge; 

else  { 

xlnitiallocal [0]  =  r[0]; 

for  (int  i  =  1;  i  <  nptles;  i++) 

xlnitiallocal [i]  =  xlnitiallocal [i-1]  +  r[i-l]  +  r[i]; 
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Appendix  C.  Practical  Extration  and  Report  Language  (PERL)  Script 

for  Parametric  Studies 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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# ! /usr/bin/env  perl 

#  This  program  automates  the  considerable  task  of  setting  up  parametric 
studies  on  the  tapered 

#  spherical  elastic  1-D  grain  problem.  It  takes  an  input  file: 
taperchain29 . cpp  and  searches 

#  through  it  replacing  the  values  of  epsilon,  q,  and  N  in  for  loops  and 
spitting  out  a  file 

#  in  the  appropriate  directory.  The  directories  are  created  on  the  fly. 

#  This  version  uses  the  updated  directories  and  is  looking  at  initial  velocity 
on  the  small 

#  grain. 

$ w  =  0.1; 

$FILE_NAME 
$SOURCE_DIR 
#$SOURCE_DIR 
$FILE_IN 
$OUT_DIR 
#$OUT_DIR 
$  E  P  S I LON_CHK 
epsilon 
$Q_CHK 
line  with  q 
$N_CHK 
with  N 

system ( "clear" ) ;  #  clear  the  screen; 

#for($w  =  0.0;  $w<=0 . 02 ;  $w+=0.01)  {  #  Restitution 

$epsilon  =1.0  -  $w;  #  taperchain . cpp  program  uses  epsilon  instead  of  w 

directly 

print "\n\n:  :::::::  w=$w  \t  epsilon  =  $epsilon  ::::::::  :\n"; 
system ("date  1 +DATE :  %m/%d/%y%nTIME :  %H:%M:%S'"); 
print " \n" ; 

chdir ( "$OUT_DIR" )  or  die  "Cant  open  $OUT_DIR" ; 
mkdir ( " w$w" ) ; 

chdir ("w$w")  or  die  "Cant  open  $OUT_DIR/w$w" ; 

for ( $N=3 ;  $N  <=20;  $N++)  { 

mkdir ( "N$N" ) ; 

chdir ("N$N")  or  die  "Cant  open  N$N" ; 

for ( $q=0 ;  $q  <=10;  $q++)  { 
mkdir ( "N$N\q$q" ) ; 

chdir ( "N$N\q$q" )  or  die  "Cant  open  N$N\q$q" ; 

$  CURRENT_D I R  =  "pwd' ; 

chop ( $CURRENT_DIR) ;  #  remove  trailing  \n 

print " Current  Dir:  $CURRENT_DIR\n" ; 

$FILE_OUT=  "$CURRENT_DIR/$FILE_NAME" ; 

open (FROM,  "$FILE_IN")  or  die  "Cant  open  $FILE_IN:  $!"; 
open (TO,  " > $ F I LE_OUT " )  or  die  "Cant  open  $CURRENT_DIR/$FILE_OUT : 

$ ! "  ; 

print" 

system ("date  1 +TIME :  %H:%M:%S'"); 

while (<FROM>)  {  #  Read  in  the  file  line  by  line  into  $_ 

#  Replacements  (Regular  expression  matching) 
if ($w  ! =  0)  { 

s/$EPSILON_CHK  \d+\ . \d+/$EPSILON_CHK  $epsilon/;  #  change 

epsilon 

} 


=  "taperchain28_w01 . cpp" ; 

=  " /home/rldoney" ; 

=  " /Users/bob/Work/Classes/Spring  2004/Dr.  Sen  study/runs"; 

=  "$SOURCE_DIR/$FILE_NAME" ; 

=  " /nf s/scratch/rldoney/TiAlV/D . SimpleTapered" ; 

=  " $SOURCE_DIR/TiAlV/D . SimpleTapered/D . Vin_large/ " ; 

=  "const  double  epsilon  =";  #  Set  pattern  to  match  line  with 

=  "const  double  q  =";  #  Set  pattern  to  match 

=  "int  nptles=";  #  Set  pattern  to  match  line 
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s/$N_CHK\d+/$N_CHK  $N/ ;  #  change  N 

s/$Q_CHK  \d+\ . \d+/$Q_CHK  $q\.0/;  #  change  q 

print  TO  $_;  #  Write  the  current  line  with  any 

changes  to  TO 

close  TO; 

} 

close  FROM; 

#  Escape  to  the  shell,  compile  the  file,  and  run  it 

#  It  was  unexpected,  but  we  need  the  1  1  because  of  the  spaces  in  the 

path 

system  ("g++  ' $FILE_OUT ' " ) ; 
system  ( " . /a . out " ) ; 
chdir ( " . . / " ) ; 

} 

chdir ( " . . / " ) ; 
print ( "\n" ) ; 

} 

system ("date  1 +DATE :  %m/%d/%y%nTIME :  %H:%M:%S'"); 
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Appendix  D.  MATLAB  Code  for  Generating  Numerical  Kinetic 

Energy  Surfaces 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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clc ;  format  compact;  format  long 
%  STC  directory 

cd  1 /Volumes/Xternal/Individual  Study/TiAlV/D . SimpleTapered/D . Vin_large 1 
pwd 

wi  =  0 ; 


qmin  =  0 ; 
qmax  =  10; 

qtot  =  qmax-qmin  +  1; 
fileprefix  =  ' taperchain_ 


for 


w=0 . 0  : 
wi  =  wi 
cd  ( [ ' w ' 
Ni  =  0; 
for  N=3 : 


.02  :  0.1 

+  1  ; 

, num2str (w) ] ) ; 


Since  N  skips,  we  need  an  index  for 
%  For  simple  tapered  chain 


N 


Build 

Build 


name 

name 


of 

of 


N#  directory 
N#q 


:  20 

Ni  =  Ni  +  1;  % 

f ilesuf f ix  =  N; 
cd  ( [ ' N ' , num2str (N) ] ) 
foldername  =  ( [ ' N ' , num2str (N) ] ) ; 
foldernameq  =  ( [foldername, ' q 1 ] ) ; 
directory 

for  q=qmin:qmax 

%qi  =  q;  %  If  we  ignore  the  monodisperse  case 

qi  =  q  +  1;  %  Indices  can't  be  zero,  so  create  qi  as  an  index 

foldername  =  ( [foldernameq, num2str (q) ]) ;  %  Foldername  is  determined 

by  current  q 

cd ( [foldername] ) 

file  =  (  [fileprefix, num2str (f ilesuf fix) ])  ; 

filename  =  (  [f ile,  1  . dat ' ] )  ;  %  Add  the  . dat  suffix  to 

filename 
pwd 


current 


%  Largest  Grain 
KEinL (Ni, qi, : ) 

E ( J)  -  KE 

%FinL (Ni , qi , : ) 
F(kN)  -  Force 


dlmread ( [filename] , ' \t ' , ' E2 . . E12 1 ) ; 
dlmread ( [filename] , ' \t 1 , ' F2 . . F1202 1 ) ; 


%  Smallest  Grain 

%t (Ni , qi , : )  =  dlmread (' taperchain_l . dat ',' \t ', 'A2 . .A5002 ') ; 

%  t (us)  -  Time 

KEoutS (Ni , qi ,  : )  =  dlmread ( ' taperchain_l . dat ' ,  ' \t '  ,  '  E2  .  . E1502 ' )  ; 

%  E(J)  -  KE 

%FoutS (Ni , qi ,  : )  =  dlmread ( 1 taperchain_l . dat 1  ,  1 \t 1  ,  1 F2  .  .  F1202  1  )  ; 

%  F(kN)  -  Force 


%  INPUTS 

o, _ 

O - 

%  Maximum  input  kinetic  energy  is  1st  element  (initial  velocity) 

%  in  the  largest  grain  of  every  chain 
KEmax_in_L (Ni , qi , : )  =  KEinL (Ni , qi , 1 ) ; 

%  Maximum  input  Force  is  somewhere  early  on  for  1st  element. 

%%  ANDed  portion  is  to  make  sure  a  small  peak  doesn't  occur  earlier 
%f or  i=l : 1200 

%  if  (FinL (Ni, qi, i+1)  <  FinL (Ni , qi , i) )  &&  (FinL (Ni, qi, i)  >  0.001) 

%  Fmax_in_L (Ni , qi )  =  FinL (Ni , qi , i ) ; 

%  break 

%  end 
%end 


%  OUTPUTS 
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%  New  technique  of  normalizing:  wrt  first  peak  hitting  last  grain. 

%  peak  occurs  when  comparing  the  difference  between  the  ith  and 
%  (i+l)st  elements.  when  it  is  less  than  zero  we  have  just  passed  a 
%  peak 

%%  ANDed  portion  is  to  make  sure  a  small  peak  doesn't  occur  earlier 
for  i=l : 1200  %  Only  need  a  small  part  of  t 

if  (KEoutS (Ni,qi, i+1)  <  KEoutS (Ni , qi , i ) )  &&  (KEoutS (Ni , qi , i )  > 

0.001) 

KEmaxF_out_S (Ni , qi )  =  KEoutS (Ni , qi , i ) ; 

break 

end 

%  For  force,  since  the  acceleration  is  the  other  way,  the  sign 
%  <  goes  to  >.  And  we  want  only  want  the  magnitude  of  the  force 
%FmaxF_out_S (Ni , qi)  =  abs (min (FoutS (Ni , qi , : ) ) ) ; 

%if  (FoutS (Ni, qi, i+1)  >  FoutS (Ni , qi , i ) ) 

%  FmaxF_out_S (Ni, qi)  =  abs (FoutS (Ni , qi , i )) ; 

%  break 
%end 
i  =  i  + 1  ; 

end 

cd  .  . 

clear  foldername 
end 

cd  .  . 

end 

%Fnorm_in_L  =  FmaxF_out_S  ./  Fmax_in_L;  %  Normalize  Fout/Fin  per 

each  specific  grain 

KEnorm_in_L  =  KEmaxF_out_S  ./  KEmax_in_L;  %  Normalize  Fout/Fin  per 

each  specific  grain 

cd  ' /' Volume s/Xternal/ Individual  Study /TiAlV/D . SimpleTapered/D . Vin_large ' 

Ni  =  linspace (3 , 20 , 18 )  ; 

qi2  =  linspace (qmin, qmax*0 . 01 , qtot) ; 

%  Figure  labeling 


if  wi  ==  1 

letter  =  '  a '  ; 
elseif  wi  ==  2 

letter  =  ' b ' ; 
elseif  wi  ==  3 

letter  =  ' c ' ; 
elseif  wi  ==  4 

letter  =  ' d ' ; 
elseif  wi  ==  5 

letter  =  ' e ' ; 

else 

letter  =  ' f  '  ; 

end 


%  For  some  reason,  these  3-d  plots  switch  the  x  and  y  plotting.  The 
%  variables  are  correct,  but  the  ordering  was  unexpected 

%plot (  squeeze (  Ft(l,2,:))  ,  squeeze (  Fin(l, 2,:)  )  )  %  N  =  3  ,  q  = 

1 

subplot (3 , 2 , wi) 

%surf (qi2 , Ni , Fnorm_in_L) 
surf (qi2 , Ni , KEnorm_in_L) 
xlabel ( 1 \bf  q 1 ) 
ylabel ( 1 \bf  N_i 1 ) 

%if  (wi  ==  1  wi  ==  3  |  wi  ==  5) 
zlabel ( ' \bf  \it  KE_N') 

%end 

text (0 , 18 , 0 . 4 , [ ' \it  \omega  =  ' , num2str (w) ] ) 


55 


%  subplot  labeling  based  on  index  of  wi 

text (0,18,0.5,  ['\bf  ( 1 , letter,  1 )  ' ]  ,  ' Font Size  1 , 14 ,  ' FontName ' ,  ' Times ' ) 

axis ([0.0  0.1  3  21  0  0.5]) 
view (135 , 10) 

caxis ( [0  0.5])  %  Adjust  the  limits  of  the  color  scheme  for  Z,  must 

come  before  colorbar 
%colorbar 

end 

cd  1 /Applications/Physics/Matlab7 1 
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Appendix  E.  Normalization 


Normalization  has  posed  some  challenges  in  trying  to  properly  assess  the  absorption  quality  of  a 
simple  tapered  chain  (STC).  Proper  normalization  schemes  will  help  the  architect  better 
determine  which  STC  is  best  for  which  application.  Recall  that  the  goal  is  to  measure  the  energy 
at  the  last  grain  versus  the  energy  put  into  the  system  by  the  first.  In  general,  the  functional  form 
will  stay  the  same  regardless  of  the  strategy  and  adjustments  in  the  normalization  will  simply 
scale  the  kinetic  energy  (KE)  surface.  In  some  cases,  the  output  force  of  each  chain  is  based  on 
the  maximum  value  felt  by  any  possible  chain  under  consideration  (i.e.,  monodisperse  and  no 
energy  loss).1  This  gives  a  measure  of  how  one  chain  is  better  than  another. 

In  this  communication,  we  have  chosen  to  form  the  ratio  based  on  the  output  KE  and  force  felt 
by  each  specific  chain.  This  serves  to  grade  the  individual  effectiveness  of  any  chain  without 
reference  to  another.  In  choosing  a  peak  value  with  this  method,  one  could  identify  either  when 
the  impulse  first  hits  the  last  grain  or  look  for  the  absolute  maximum  peak  whenever  it  may 
occur.  We  have  chosen  to  use  the  former  for  several  reasons.  First,  it  ignores  the  complexity  of 
nonlinear  reverberations  which  can  lead  to  large  peaks  at  unpredictable  times.  Second,  we  argue 
that  this  is  just  as  realistic  as  selecting  the  maximum  value  anywhere  in  the  time  spectrum. 

It  turns  out  that  in  most  cases,  the  absolute  maximum  is  the  first  peak.  There  are  special  cases 
where  the  maximum  may  occur  at  later  times  and  this  needs  to  be  investigated  further.  In 
figure  6  (in  the  body  of  the  report)  for  q  =  0.1,  for  example,  we  see  the  striking  occurrence  of  the 
secondary  pulse  about  225  ps  being  much  stronger  than  the  initial  arrival  at  35  ps.  This  is  one  of 
those  instances  that  disagrees  with  the  way  we  choose  to  normalize  our  KE  surfaces.  We  have 
investigated  this  particular  case  further  without  including  extraneous  plots  and  report  the 
following  observations.  The  effect  exists  for  N  =  15  -20  for  constant  co.  When  N  =  20  is  held 
fixed  and  restitution  is  increased,  the  peak  KE  once  again  occurs  for  the  first  arrival  of  the  pulse. 
As  q  increases,  so  do  the  number  of  collisions  and  the  requisite  energy  loss  (since  co  ^  0). 
Therefore,  it  is  less  likely  to  find  a  global  peak  later  in  the  simulation.  The  situation  is  further 
complicated  by  the  interplay  between  q  and  co. 
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ADELPHI  MD  20783-1197 


ABERDEEN  PROVING  GROUND 


1  DIR  USARL 

AMSRD  ARL  Cl  OK  TP  (BLDG  4600) 
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NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


5  INST  FOR  DEFENSE  ANALYSES 
H  BERTRAND 
G  BOEZER 
J  HEAGY 
I  KOHLBERG 
M  RIGDON 

4850  MARK  CENTER  DR 
ALEXANDRIA  VA  22311-1882 

7  STATE  UNIVERSITY  OF  NEW  YORK 
DEPT  OF  PHYSICS 
R  GONSALVES 
B  POWELL 
S  SEN  (5  CPS) 

239  FRONCZIAK  HALL 
BUFFALO  NY  14260-1500 


AMSRD  ARL  WM  TC 
M  FERMEN-COKER 
AMSRD  ARL  WM  TE 
P  BERNING 
C  HUMMER 
J  POWELL 
AMSRD  ARL  SL  BE 
B  BRUCHEY 
T  BJERKE 


1  STATE  UNIVERSITY  OF  NEW  YORK 
DEPT  OF  MATHEMATICS 
B  PITMAN 

244  MATHEMATICS  BLDG 
BUFFALO  NY  14260-2900 


ABERDEEN  PROVING  GROUND 


37  DIR  US  ARL 

AMSRD  ARL  WM 
T  WRIGHT 

AMSRD  ARL  WM  MD 
B  CHEESEMAN 
AMSRD  ARL  WM  T 
B  BURNS 
M  ZOLTOSKI 
AMSRD  ARL  WM  TA 
A  BARD 

P  BARTKOWSKI 
R  DONEY  (10  CPS) 
G  FILBEY 
B  FREY 

V  HERNANDEZ 
E  HORWATH 
T  JONES 
M  KEELE 
P  KINGMAN 
D  KLEPONIS 
C  KRAUTHAUSER 
B  LEAVY 
I  LISKO 
A  MIHALCIN 
J  RUNYEON 
S  SCHOENFELD 
AMSRD  ARL  WM  TB 
R  SKAGGS 
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Intentionally  left  blank. 
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